micrOMEGAs and the relic density in the MSSM 



G. Belanger 

LAPTH,C/iemm de Bellevue, B.P. 110, F-74941 Annecy-le-Vieux, Cedex, France. 

Abstract 

micrOMEGAs is a program that calculates the relic density of the lightest supersym- 
metric particle (LSP) in the minimal super symmetric standard model. All tree-level 
processes for the annihilation of the LSP are included as well as all possible coanni- 
hilation processes. The cross-sections extracted from CompHEP are calculated exactly. 
Relativistic formulae for the thermal average are used and care is taken to handle 
poles and thresholds by adopting specific integration routines. The Higgs masses are 
calculated with FeynHiggsFast and the QCD corrected Higgs widths with HDECAY. 



1 Introduction 

One of the strong arguments in favour of supersymmetry is that R-parity conserving super- 
symmetric models have a cold dark matter candidate(CDM), the lightest supersymmetric 
particle(LSP). The preferred candidate is a neutralino. The contribution of the LSP to the 
relic density is however very model dependent and varies by several orders of magnitude over 
the whole allowed parameter space of the mininal supersymmetric standard model (MSSM). 
The measurements of the relic density then imposes stringent constraints on the parameters 
of the MSSM often favouring solutions with light supersymmetric particles. There are basi- 
cally two mechanisms that can significantly reduce the estimate of the relic density leading 
to acceptable values while having a not so light supersymmetric spectrum: annihilation 
via a s-channel resonance and coannihilations where the LSP interacts with slightly heavier 
sparticles. Special care must be taken to treat these two cases carefully if one wants to make 
predictions for the relic density of CDM at the few percent level. Although, at present, the 
generally assumed range for the relic density of CDM is .1 < Qh"^ < .3 , the recent data 
from BOOMERANG [Ul coupled with some constraints already indicates a more restrictive 
range at 2a, .11 < flh"^ < .150. As new measurements will be performed in the near future, 
improvements over the present limits are expected. 

There exist many calculations of the relic density in the MSSM using various approx- 
imations both for the evaluation of the thermally averaged cross-section and for solving 
the Boltzmann equation for the density of dark matter particles [3-8]. Among these, 
Neutdriver0] and DarkSusy[§] are publicly available. Our purpose was to provide a tool 
that evaluates with high accuracy the annihilation cross-sections even in regions near poles 
and thresholds, that is both flexible and upgradable and that goes beyond DarkSusy as far 
as the calculation of the relic density is concerned ||]. 

The first calculations of the relic density used an expansion of the annihilation cross- 
section in power series of neutralino velocities, a = a + bv"^. While this approximation 
works well in large regions of parameter space, it fails when annihilation through s-channel 
resonance is important. When the masses are such that the neutralinos can annihilate 
through a s-channel Higgs resonance |^ |TU[ or a s-channel Z resonance^, the annihilation 



rate increases significantly often bringing the relic density in an acceptable range as one gets 



close to the s-channel pole. In fact as one gets very close to the pole, the annihilation rate 
becomes so fast that the neutralinos cannot constitute the only source of dark matter. These 
effect are especially important in models with large tan (3 where the heavy Higgs resonances 
are very wide. The proper relativistic formalism for the evaluation of the thermally averaged 
cross-section was then introduced by |TI|. The relativistic formalism, generalized to the case 

was implemented in the codes of and 



of coannihilations iT2 



for the case of gaugino 

coannihilations. We follow basically this formalism, although contrary to DarkSusy, we still 
rely on approximations for the solution of the relic density equations and the determination 
of the freeze-out temperature. This allows to significantly increase the speed of the program 
and proves to be very useful when scanning over a large parameter space. 

Coannihilation processes where the LSP interacts with slightly heavier sparticles can oc- 
cur in principle with any supersymmetric particle [|13|, although in SUGRA models, the most 
common coannihilations are with gauginos or right-handed sleptons. The importance of the 



coannihilation channels were emphasized before both for gauginos |T^, |T^, sleptons [|T5], |T6 
or stops |18|. In micrOMEGAs we include ALL coannihilation channels, in all more than 
2800 processes not counting charged conjugate processes. The tree-level cross-sections are 
calculated exactly including the full set of diagrams contributing to each process. The calcu- 
lations of cross-sections are based on CompHEP|]T^, an automatic program for the evaluation 
of tree-level Feynman diagrams. Furthermore we include also some higher order effects. 



namely the two- loop corrections to the Higgs masses |]20[ and the one- loop corrections to 
the Higgs widths[^]. The latter turns out to be important in the large tanjS region. Al- 
though we will generally assume that the neutralino is the LSP, micrOMEGAs can be used 
to compute the relic density with any supersymmetric particle as the LSP, in particular the 
sneutrino. This is because all (co-) annihilation of any pairs of supersymmetric particles into 
any pairs of standard model or Higgs particles are included. 

The main characteristics of micrOMEGAs, are 

• Complete tree-level matrix elements for all subprocesses 

• Includes all coannihilation channels with gauginos, sleptons and squarks. 

• Loop-corrected Higgs masses and widths 

• Speed of calculation 

After the important equations for the calculation of the relic density are summarized, we 
give a short description of the parameters of the supersymmetric model and of the package. 
Finally we present some results and comparisons with another program in the public domain, 
DarkSusy. 



2 Calculation of the relic density 

The relic density at present is calculated from 

n^oh^ = 2.755 X 10^-7^7^0 (2.1) 

^1 GeV 

where Yq is the abundance of the LSP today. To find Yq = Y{T = Tq), one needs to solve 
the evolution equation for Y 



dY_ 7cg,{T) 
dT~\ AhG 



<av> {Y^ - Yl) (2.2) 



gtf{T) is a degrees of freedom parameter derived from the thermodynamics describing the 



state of the universe ^ and Y^q = Yeq{T) represents the thermal equihbrium abundance 



where we sum over all supersymmetric particles i with mass rrij and Qi degrees of freedom. 
Kn is the modified Bessel function of the second kind of order n . Note that Yeq falls 
rather rapidly as the temperature decreases. < af > is the relativistic thermally averaged 
annihilation cross-section 

E9i9j I ds^Ki{^/T)p1-a,j{s) 
< ov >= ^"^'-^'"-^^ ^ , , (2.4) 

i 

where aij is the total cross section for annihilation of a pair of supersymmetric particles 
into some Standard Model particles, and Pij is the momentum of the incoming particles in 
their center-of-mass frame. The summation is over all supersymmetric particles. Integrating 
Eq. p]^ from T = oo to T = Tq would lead Yq. 



Although one can solve for Y numerically, the procedure is extremely time consuming 
especially when scanning over a large parameter space and when we include a great number 
of processes. It is therefore important to seek as good an approximation as possible to 
speed up the code. We follow the usual procedure of defining a freeze-out temperature 
T/[0. which can be extracted by solving iteratively 



dln{Y,q) hgJJl ^ v Aa ^ 9^ ^9 

where 6 is some small constant number. The freeze-out temperature is defined from Yf = 
Y{Tf) = (1 + 6)Yeq(Tf). Starting from a typical value Tf = only a few iterations 

are necessary to find a solution to this equation with 5 = 1.5 ± .2. 

In the second regime, where Y ^ Yeq, one can neglect Y^^ completely. One finds [pi]] 

Furthermore we find that the solution ( ^^.6] ) does not depend significantly on 6. 

Typical freeze-out temperatures vary between IGeV and lOGeV, in this temperature 



range, the h^ff and functions can vary by about 20% [0. To achieve an accuracy better 
than 10%, one cannot use a constant value for these functions. We use the numerical tables 
of the DarkSUSY package 0] for a precise evaluation of heff and g^:. 



2.1 Numerical integration and summation 

In order to find Tf by solving Eq. |2.5| , we have to evaluate several integrals and perform a 
summation over different annihilation channels. In the evaluation of the thermally averaged 
cross-section we have included all two-body subprocesses involving two LSP's, the LSP and 
a co-annihilating SUSY particle (all the particles of the MSSM) as well as all subprocesses 



involving two coannihilating SUSY particles. The final states include all possible standard 
model and Higgs particles that contribute to a given process at tree-level. The total number 
of processes exceeds 2800 not including charged conjugate processes. In practice processes 
involving the heavier SUSY particles contribute only when there is a near mass degeneracy 
with the LSP since there is a strong Boltzmann suppression factor -B/ in Eq. |2.4| . 

where mi,mj are the masses of the incoming particles. To speed up the program a given 
subprocess is removed from the sum ( [2.4] ) if the total mass of the incoming particles is such 
that Bf is below some limit, defined by the user. Alhtough it would be sufficient to 
use B^ = 10~^ to give a precision of 1% when acoan ~ '^x?x?' ^ more restrictive 

value, B^ = 10'^, to allow for cases where acoan ^ '^x?x?' '^^^^ occur for example for 
coannihilation processes with squarks which depend on as, for processes with poles or in 
some regions of parameter space where cr^o^o is suppressed. 

In our program we provide two options to do the integrations, the fast one and the 
accurate one. The fast mode already gives a precision of about 1% which is good enough for 
all practical purposes. In the accurate mode the program evaluates all integrals by means of 
an adaptative Simpson program. It automatically detects all singularities of the integrands 
and checks the precision. In the case of the fast mode the accuracy is not checked. We 
integrate the squared matrix elements over cos 6', the scattering angle, by means of a 5 



point Gauss formula. For integration over s, Eq. |2]J we use a restricted set of points which 
depends whether we are in the vicinity of a s-channel Higgs(Z,W) resonance or not. We 
increase the number of poins if the Boltzmann factor corresponding to rripoie is larger than 
0.01 -R. 



3 MSSM parameters, Higgs mass and widths 

The input parameters are the ones of the soft SUSY Lagrangian defined at the weak scale. 



using the same notation as in CompHEP/SUSY models[Q. In the model used, the masses 
of fermions of the first generation are set to zero. The masses of the quarks of the second 
generation are also set to zero, so that there is no mixing of the squarks of the first two 
generations. However we have kept the mass for the muon as well as the trilinear coupling 
as this is relevant for the calculation of the muon anomalous magnetic moment. Although 
first generation sleptons are pure left /right states, they are properly ordered according to 
their masses so that the correct coannihilating particle corresponding to the lightest slepton 
is taken into account. The sign convention for the parameters fi and A, is given explicitly 
in i. 

The calculation of the Higgs masses are done with FeynHiggsFast |20[] . For the Higgs 



widths, it is necessary to take into account the QCD corrections to the partial widths, 
h{H, A) bb. This is particularly important at large tan (3 where the bb mode is the 
dominant one. The QCD corrections can be very large for heavy Higgses, easily a factor 
of two above the tree-level value for a Higgs of ITeV, due mostly to the running of the 
quark mass at high scale. To take these corrections into account we have redefined the 
vertices hqq, Hqq and Aqq using an effective mass that reproduces the radiatively corrected 



width. We used HDECAY|21[] to produce a table of mass-dependent QCD-corrected Higgs 
partial widths. From this, the effective quark masses mb(mj^) are extracted and simple 
interpolation can reproduce the one- loop corrected width for any value of mn- In the region 
of physical interest, the precision on the width is at the per- mil level safe for the neutral 
Higgses partial widths near the ti threshold. However, this region does not contribute 
significantly to the neutralino cross-section. For the charged Higgs, one can extract from 
HDECAY both an effective rrit as well as an effective rrib using high and low tan/] values. 
This way we could reproduce at better than 1% level the partial width — > tb. The 
effective quark mass mb{Q) evaluated a.t Q = mi + 1712, the sum of the initial masses, is 
used by default in the {h, H, A) — bb vertices. This will lead to the correct result when the 
contribution of the Higgs resonance is very important {Mh ~ 2m^o). 



4 Description of micrOMEGAs 

micrOMEGAs is a C program that also calls some external FORTRAN functions. micrOMEGAs 
relies on CompHEP[|T^ for the definition of the parameters and the evaluation of all cross- 
sections. The program is contained in a package that lets the user choose between weak scale 
parameters or parameters of SUGRA models as input parameters. The latter is achieved 
through a link with I SASUGRA/ 1 s a ] e t [^ ] . The package can be obtained at 



http : / /wwwlapp . in2p3 . f r/lapth/niicromegas| 



As we have already mentionned, due to a strong Boltzmann suppression factor, only a 
small fraction of the available processes are needed, those with a sparticle close in mass to the 
LSP. In principle, compilation of the full set of subprocesses is possible, but such a program 
would be huge and could not be distributed easily. To avoid this problem, we include in our 
package the program CompHEP[|l^] which generates, while running, the subprocesses needed 
for a given set of MSSM parameters. The generated code is linked during the run to the 
main program and executed. The corresponding "shared" library is stored on the user disk 
space and is accessible for all subsequent calls, thus each process is generated and compiled 
only once. Such approach can be realized only on Unix platforms which support dynamic 
linking. 

The complete list of input soft SUSY parameters at the weak scale can be found in 
0. When using the SUGRA option, the 7 input parameters are the usual 5 parameters of 
a SUGRA model, defined at the GUT scale, mo,mi/2,Ao,ta.nP,sign{fi). In addition one 
must specify the value of the top mass as well as the the Isajet option (model = 1, 2 ) for 
a SUGRA model with or without gauge coupling unification at the GUT scale. When this 
option is used, the value of the weak scale soft supersymmetric Lagrangian are then extracted 
from ISASUGRA. Whether the input parameters are defined at the weak scale or at the GUT 
scale, one can always choose to redefine additional parameters, such as the standard model 
parameters or the widths of Higgses or SUSY particles. In processes with t-channel poles it 
is sometimes necessary to specify a width for some particles such as gauginos. By default 
these widths have been set to IGeV. Even though the LSP is assumed to be stable, to avoid 
any spurious pole it is necessary to introduce also a small width. Numerically, the results 
do not depend on the exact value chosen for these widths. 

After reading the input parameters and calculating the physical parameters needed for 
the evaluation of the cross-sections using the functions of CompHEP, the calculation of the 
relic density is performed. The value of flh^ as well as the list of channels that give the 



most significant contribution to it are given. For all specifications on the functions available 
and the options that can be set see 0. In addition we provide subroutines that calculate 
various constraints on the MSSM parameters: direct limits from colliders, Ap, b s'j and 
{g — 2)^. All these constraints can be updated or replaced easily. 

5 Results and Comparisons 

Table 1: Sample results and comparison with DarkSusy 
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The micrOMEGAs code was extensively tested against another public package for calcu- 
lating the relic density, DarkSUSY. As discussed previously, the two codes differ somewhat 
in the numerical method used for solving the density equations. micrOMEGAs includes more 
subprocesses(e.g. all coannihilations with sfermions), loop-corrected Higgs widths, and com- 
plete tree-level matrix elements for all processes. DarkSUSY includes on the other hand some 



loop induced processes such as XX ~^ 99^ 77 which are generally small. Whenever the coan- 
nihilation channels with sfermions and the Higgs pole are not important we expect good 
agreement with DarkSUSY. We have first compared micrOMEGAs with a version of DarkSUSY 
where we have replaced the matrix elements by CompHEP matrix elements. As expected, 
we found excellent numerical agreement between the two programs (at the 2% level) for all 
points tested. We have then made comparisons with the original version of DarkSUSY. The 
results of these comparisons are displayed for a few test points in Table |l]. We found in gen- 
eral good agreement between the two codes. However we have observed some discrepancies 
that could reach up to (30%) in particular in the process XiXi ~^ t5(model B in Table 0)|5 
As displayed in the line ^^reei when removing non-gaugino coannihilation channels and re- 
verting to the tree-level treatment of the Higgs width we recover results similar to DarkSUSY. 
The impact of these extra channels, model C for sleptons and model D for squarks can be 
as large as an order of magnitude and depends critically on the mass difference with the 
lightest neutralino 0. 



Figure 1: Vlh? vs the NLSP-LSP mass difference for a) model C with Mli(Mri) as a free 
parameter. The NLSP is a fi (z/); b) model D with Mui(Mdi) as a free parameter. The 
NLSP is the ti {hi). For all cases, the value of Vth"^ neglecting coannihilation channels is 
also shown (dotted lines). 
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In Fig. |l] the variation of the relic density as function of the mass difference between 
the neutralino x? and the Next-to-Lightest Supersymmetric particle(NLSP) is displayed. In 
Fig. the SUSY parameters correspond to model C, with either the right-handed slepton 
masses Mrl=Mr2=Mr3 or the left-handed slepton masses M11=M12=M13 as free parameters 
(in this case we have also fixed Mr3=289 GeV). For the former choice of parameters, the 
f is the NLSP while for the latter the NLSP is a sneutrino. Model C corresponds to 
AMf_^ = 10.4GeV. In a model where the neutralino is the only source of darkmatter, the 
lower bound on the relic density then implies a minimum mass difference between the LSP 
and the NLSP(> 5 — lOGeV), avoiding the difficulties of detecting a nearly degenerate NLSP 
at colliders. In Fig.|I]b, the dependence of the relic density on the mass difference between 



-'^The CompHEP result for this matrix element agrees with the result of GraceSUSY|26 



^Extensive comparisons of micrOMEGAs with an improved version of DarkSusy including slepton coan- 
nihilation channels were also performed recently in p8[, complete agreement was found. 



a coanihilating squark and the LSP is displayed. The SUSY parameters are those of model 
D with either the right-handed u-squark or d-squark masses kept as free parameters. These 
two cases lead to a t or a 6 NLSP respectively. Model D corresponds to AMi_^ = 26.9GeV. 
Due to the large cross-sections in coannihilation channels involving strongly interacting 
particles, One witnesses a sharp drop in the relic density as soon as AM^lsp-lsp drops 
below 50GeV. 



Fi^^e 2: Comparison of ilh^ calculated liS^g 
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As stressed above, the effect of the Higgs width is particularly important at large tan f3 
with the enhanced contribution of the 6-quark coupling to the heavy scalar Higgs. However, 
the one-loop QCD corrections to the widths amount to a reduced effective b-quark mass 
and a much smaller width especially at large values of rriH- If it was not for the strong 
Boltzmann suppression factor singling out the contribution at ^/s ~ 2m^ there will be little 
difference after integrating over the peak for the one- loop or tree- level result. However the 
effect observed can be as much as a factor 2. For ~ Ma/ '2^ the narrower resonance (1- 
loop) suffers less from the Boltzmann reduction factor leading to < a^-^°°^ > / < <^tree >> 1 
and VLi-ioop < ^tree- Further away from the pole however one catches the contribution from 
the wider resonance without excessive damping from the Boltzmann factor, as expected 
^i-ioop > ^tree- This is illustrated both for the parameters of models E and F while varying 
Ma in a wide range around the values Ma = 2m^. 

Our numerical results were also compared with Ref. |p7||. Qualitative agreement is found 
in the case of a SUGRA model although we use a different RGE code. Precise comparisons 
necessitates a careful tuning of parameters to make sure we have the same parameters at 
the weak scale. A random scan over mo — mi/2 for /i > shows the typical shape of the 
allowed region (.1 < Qh'^ < .3) in SUGRA models for moderate values of tan/3. For large 
tan i3, one clearly see the dramatic effect of the heavy Higgs pole in the central band that 
goes up to very large values of mo, m,i/2- Anywhere in this band or in the coannihilation 
tails a generally heavy supersymmetric particles spectrum can be expected. 



Figure 3: Qh"^ in the mQ—mi/2 plane in a SUGRA model with = 0, /i > 0, rritop = 175GeV 
and a)tan/3 = 10 b)tan/? = 50. The dark region corresponds to .1 < Qh"^ < .3 and the light 
(green) region to Qh"^ < .1. In the white region at large mo — Tni/2, the relic density is above 
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6 Conclusion 

The package micrOMEGAs allows to calculate the relic density of the LSP in the most general 
MSSM with Rp conservation. This is the first program that includes all possible coannihila- 
tion channels^. The package is self-contained safe for the ISASUGRA package that is required 
when using the SUGRA option. All possible channels for coannihilations are included and 
all matrix elements are calculated exactly at tree level with the help of CompHEP. Loop 
corrections for the masses of Higgs particles (two-loop) and the width of the Higgs (QCD 
one-loop) are implemented. Good agreement with existing calculations is found when iden- 
tical set of channels are included. Future versions will include interfaces to other codes that 
use the renormalization group equations to calculate the weak scale parameters. Including 
loop corrections to neutralino masses is also planned. Even though these corrections are 
only a few GeV's they can alter significantly the calculation of the relic density when there 
is a near mass degeneracy with the next to lightest supersymmetric particle that contributes 



to a coannihilation channel Although the loop processes are in general small, we plan 



to include XX ~^ 77) 99 iii cin update of micrOMEGAs. 
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